Spatial extent of an outbreak in animal epidemics 
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Characterizing the spatial extent of epidemics at the outbreak stage 
is key to controlling the evolution of the disease. At the outbreak, 
the number of infected individuals is typically small, so that fluctu- 
ations around their average are important: then, it is commonly as- 
sumed that the susceptible-infected-recovered (SIR) mechanism can 
be described by a stochastic birth-death process of Galton-Watson 
type. The displacements of the infected individuals can be modelled 
by resorting to Brownian motion, which is applicable when long- 
range movements and complex network interactions can be safely 
neglected, as in case of animal epidemics. In this context, the spa- 
tial extent of an epidemic can be assessed by computing the convex 
hull enclosing the infected individuals at a given time. We derive the 
exact evolution equations for the mean perimeter and the mean area 
of the convex hull, and compare them with Monte Carlo simulations. 
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Models of epidemics traditionally consider three classes of 
populations, namely, the susceptibles (S), the infected 
(I), and the recovered (R). This provides the basis of the so- 
called SIR model [T] [3], a fully connected mean-field model 
where the population sizes of the three species evolve with 
time t via the coupled nonlinear equations: dS/dt = ~j3IS; 
dl jdt = f3IS—jI and dR/dt = jl. Here 7 is the rate at which 
an infected individual recovers and /3 denotes the rate at which 
it transmits the disease to a susceptible [4j [5] E] ■ In the sim- 
plest version of these models, the recovered can not be infected 
again. These rate equations conserve the total population size 
I(t) + S(t) + R(t) — N, and one assumes that initially there 
is only one infected individual and the rest of the population 
is susceptible: 7(0) = 1, S(0) = N - 1, and R(0) = 0. Of 
particular interest is the outbreak stage, i.e., the early times 
of the epidemic process, when the susceptible population is 
much larger than the number of infected or recovered. Dur- 
ing this regime, for large N, the susceptible population hardly 
evolves and stays S(t) ~ N, so that nonlinear effects can be 
safely neglected and one can just monitor the evolution of the 
infected population alone: dl/dt~ (PN — ~/)I(t). Thus, the ul- 
timate fate of the epidemics depends on the key dimensionless 
parameter Ro = /3./V/7, which is called the reproduction rate. 
If Rq > 1 the epidemic explodes and invades a finite fraction 
of the population, if Ro < 1 the epidemic goes to extinction, 
and in the critical case Ro — 1 the infected population remains 

constant [7][8l|9]. 

This basic deterministic SIR has been generalized to a va- 
riety of both deterministic as well as stochastic models, whose 
distinct advantages and shortcomings arc discussed at length 
in [3] 1101 [TT] . Generally speaking, stochastic models are more 
suitable in presence of a small number of infected individuals, 
when fluctuations around the average may be relevant [3] [10] . 
During the outbreak of epidemics, the infected population is 
typically small: in this regime, the evolution can be modeled 
by resorting to a stochastic birth-death branching process of 
the Galton-Watson type for the number of infected [3l llO|[TT| . 
where each infected individual transmits the disease to an- 
other individual at rate /3N and recovers at rate 7. The epi- 
demic may become endemic for Ro > 1, becomes extinct for 
Ro < 1, whereas for Ro = 1 fluctuations are typically long 
lived and completely control the time evolution of the infected 
population |4j [5] [6] . 



How far in space can an epidemic spread? Branching pro- 
cesses alone are not sufficient to describe an outbreak, and 
spatial effects must necessarily be considered [T1I51 I12II13II14] . 
Quantifying the geographical spread of an epidemic is closely 
related to the modelling of the population displacements. 
Brownian motion is often considered as a paradigm for de- 
scribing the migration of individuals, despite some well-known 
shortcomings: for instance, finite speed effects and preferen- 
tial displacements are neglected. Most importantly, a num- 
ber of recent studies have clearly shown that individuals ge- 
ographically far apart can actually be closely related to each 
other through the so-called small-world connections, such as 
air traffic, public transportation and so on: then, the spread 
of epidemics among humans can not be realistically modelled 
without considering these complex networks of interconnec- 
tions [151 1161 im I18| . Nonetheless, Brownian motion provides 
a reasonable basis for studying disease propagation in animals 
and possibly in plants (here, pathogen vectors are insects) f4]. 

While theoretical models based on branching Brownian 
motion have provided important insights on how the pop- 
ulation size grows and fluctuates with time in a given do- 
main [Tl 151 IT31 114] . another fundamental question is how the 
spatial extension of the infected population evolves with time. 
Assessing the geographical area travelled by a disease is key 
to the control of epidemics, and this is especially true at the 
outbreak, when confinement and vaccination could be most 
effective. One major challenge in this very practical field of 
disease control is how to quantify the area that needs to be 
quarantined during the outbreak. For animal epidemics, this 
issue has been investigated experimentally, for instance in the 
case of equine influenza [19]. The most popular and widely 
used method for this consists in recording the set of positions 
of the infected animals and, at each time instant, construct a 
convex hull, i. e.,a minimum convex polygon surrounding the 
positions (Fig. [I] for a precise definition of the convex hull, 
see below). The convex hull at time t then provides a rough 
measure of the area over which the infections have spread up 
to time t. The convex hull method is also used to estimate 
the home range of animals, i.e., the territory explored by a 
herd of animals during their daily search for food [201 121| . 

In this paper, we model the outbreak of an epidemic as 
a Galton-Watson branching process in presence of Brownian 
spatial diffusion. Despite infection dynamics being relatively 
simple, the corresponding convex hull is a rather complex 
function of the trajectories of the infected individuals up to 
time t, whose statistical properties seem to be a formidable 
problem. Our main goal is to characterize the time evolution 
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of the convex hull associated to this process, in particular its 
mean perimeter and area. 

The rest of the paper is organized as follows. We first de- 
scribe precisely the model and summarize our main results. 
Then, we provide a derivation of our analytical findings, sup- 
ported by extensive numerical simulations. We conclude with 
perspectives and discussions. Some details of the computa- 
tions are relegated to the Supplementary Material. 

The model and the main results 

Consider a population of N individuals, uniformly distributed 
in a two dimensional plane, with a single infected at the ori- 
gin at the initial time. At the outbreak, it is sufficient to keep 
track of the positions of the infected, which will be marked as 
'particles'. The dynamics of the infected individuals is gov- 
erned by the following stochastic rules. In a small time interval 
dt, each infected alternatively 

(i) recovers with probability 'ydt. This corresponds to the 
death of a particle with rate 7. 

(ii) infects, via local contact, a new susceptible individual from 
the background with probability bdt. This corresponds to the 
birth of a new particle that can subsequently diffuse. The 
originally infected particle still remains infected, which means 
that the trajectory of the originally infected particle branches 
into two new trajectories. The rate b replaces the rate (3N in 
the SIR or the Galton- Watson process mentioned before. 

(iii) diffuses with diffusion constant D with probability I — 
(7 + b) dt. The coordinates {x(t),y(t)} of the particle get 
updated to the new values {x(t) + r/ x (t) dt,y(t) + rj v (t)dt}, 
where rj x (t) and r] y (t) are independent Gaussian white noises 
with zero mean and correlators {r] x (t)T]x(t')) = 2D8(t — t'), 
(%(*)%(*')> = 2D8(t - t') and (*)%(*')} = 0. 

The only dimensionless parameter in the model is the ratio 
Ro = 6/t> i- e > t ne basic reproduction number. 

Consider now a particular history of the assembly of the 
trajectories of all the infected individuals up to time t, start- 
ing from a single infected initially at the origin (see Fig. [I]). 
For every realization of the process, we construct the asso- 
ciated convex hull C. To visualize the convex hull, imagine 
stretching a rubber band so that it includes all the points of 
the set at time t inside it and then releasing the rubber band. 
It shrinks and finally gets stuck when it touches some points 
of the set, so that it can not shrink any further. This final 
shape is precisely the convex hull associated to this set. 

In this paper, we show that the mean perimeter (L(t)) and 
the mean area (A(t)) of the convex hull are ruled by two cou- 
pled nonlinear partial differential equations that can be solved 
numerically for all t (see Fig. [2|. The asymptotic behavior for 
large t can be determined analytically for the critical (Ro = 1), 
subcritical (Ro < 1) and supercritical (Ro > 1) regimes. In 
particular, in the critical regime the mean perimeter saturates 
to a finite value as t — > 00, while the mean area diverges log- 
arithmically for large t 
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This prediction seems rather paradoxical at a first glance. 
How can the perimeter of a polygon be finite while its area 
is divergent? The resolution to this paradox owes its origin 
precisely to statistical fluctuations. The results in Eqs. I] 
and [2] are true only on average. Of course, for each sample, 
the convex hull has a finite perimeter and a finite area. How- 



ever, as we later show, the probability distributions of these 
random variables have power-law tails at long time limits. For 
instance, while Prob(L, t — > 00) ~ L~ 3 for large L (thus lead- 
ing to a finite first moment), the area distribution behaves as 
Prob(A, t — > 00) ~ A~ 2 for large A. Hence the mean area is 
divergent as t — > 00 (see Fig. [2]). 

When Ro =fc 1, the evolution of the epidemic is controlled 
by a characteristic time t* , which scales like t* ~ |_Ro — 
For times t <t* the epidemic behaves as in the critical regime. 
In the subcritical regime, for t > t* the quantities (L(t)) and 
(A(t)) rapidly saturate and the epidemic goes eventually to ex- 
tinction. In contrast, in the supercritical regime (which is the 
most relevant for virulent epidemics that spread fast), a new 
time-dependent behavior emerges when t > t* , since there ex- 
ists a finite probability (namely 1 — l/Ro) that epidemic never 
goes to extinction (Fig. pj. More precisely, we later show that 

(L(t>t*)) = 4n{l-±-) y/Dj(Rc-l)t [3] 



(A(t»f)> = 4tt (l-^J D-y(R 



[4] 



The ballistic growth of the convex hull stems from an un- 
derlying travelling front solution of the non-linear equation 
governing the convex hull behavior. Indeed, the prefactor 
of the perimeter growth is proportional to the front veloc- 
ity v* — 2\J D 7 (Ro — 1). As time increases, the susceptible 
population decreases due to the growth of the infected individ- 
uals: this depletion effect leads to a breakdown of the outbreak 
regime and to a slowing down of the epidemic propagation. 

The statistics of the convex hull 

Characterizing the fluctuating geometry of C is a formidable 
task even in absence of branching (b = 0) and death (7 = 0), 
i.e., purely for diffusion process in two dimensions. Major 
recent breakthroughs have nonetheless been obtained for dif- 
fusion processes [23H24| by a clever adaptation of the Cauchy's 
integral geometric formulae for the perimeter and area of any 
closed convex curve in two dimensions. In fact, the prob- 
lem of computing the mean perimeter and area of the convex 
hull of any generic two dimensional stochastic process can be 
mapped, using Cauchy's formulae, to the problem of comput- 
ing the moments of the maximum and the time at which the 
maximum occurs for the associated one dimensional compo- 
nent stochastic process [23 , 24 . This was used for computing, 
e.g., the mean perimeter and area of the convex hull of a two 
dimensional regular Brownian motion [231 124] and of a two 
dimensional random acceleration process |25| . 

Our main idea here is to extend this method to compute 
the convex hull statistics for the two dimensional branching 
Brownian motion. Following this general mapping and using 
isotropy in space (see Supplementary Materials), the average 
perimeter and area of the convex hull are given by 



(L(t)) 
(A(t)) 



2n{x m (t)) 



(v(t m ))] 



where x m is the maximum displacement of our two- 
dimensional stochastic process in the x direction up to time 
t, t m is the time at which the maximum displacement along 
x direction occurs and y(t m ) is the ordinate of the process at 
t m , i.e., when the displacement along the x direction is maxi- 
mal. A schematic representation is provided in Fig. |4j where 
the global maximum x m is achieved by one single infected in- 
dividual, whose path is marked in red. A crucial observation 
is that the y component of the trajectory connecting O to 
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this red path is a regular one dimensional Brownian motion. 
Hence, given t m and t, clearly {y 2 (tm)} = 2D(t m )- Therefore, 



<A(t)> = 7T \(x 2 m (t)) - 2D{t m (t))} 



[7] 



Equations [5] and thus show that the mean perimeter and 
area of the epidemics outbreak are related to the extreme 
statistics of a one dimensional branching Brownian motion 
with death. Indeed, if we can compute the joint distribution 
Pt(x m ,tm), we can in turn compute the three moments (x m ), 
(x'm} and (t m ) that are needed in Eqs. [H] and [7]. We show 
below that this can be performed exactly. 

The convex hull perimeter and the maximum X m .For the 
average perimeter, we just need the first moment (x m (t)) = 
J °° x m qt(x m ) dx m , where qt(x m ) denotes the probability den- 
sity of the of the maximum of the one dimensional com- 
ponent process. It is convenient to consider the cumula- 
tive distribution Qt(x m ) i.e., the probability that the max- 
imum ^-displacement stays below a given value x m up to 
time t. Then, q t (x m ) = dQ t (x m )/dx m and (x m (t)) = 
Jn°°[l ~~ Qt(x m )] dx m . Since the process starts at the origin, 
its maximum a;- displacement, for any time t, is necessarily 
nonnegative, i.e., x m > 0. We next write down a backward 
Fokker-Planck equation describing the evolution of Qt(x m ) by 
considering the three mutually exclusive stochastic moves in 
a small time interval dt: starting at the origin at t = 0, the 
walker during the subsequent interval [0, dt] dies with prob- 
ability ydt, infects another individual (i.e., branches) with 
probability bdt = Ro'ydt, or diffuses by a random displace- 
ment Ax — r] x (0) dt with probability 1 — 7(1 + Ro)dt. In the 
last case, its new starting position is Ax for the subsequent 
evolution. Hence, for all x m > 0, one can write 

Qt+dt(x m ) = -ydt + R ^dtQ 2 t (x m ) 
+[l-j(Ro + l)]dt(Q t (x m -Ax)), [8] 

where the expectation (} is taken with respect to the random 
displacements Aa;. The first term means that if the process 
dies right at the start, its maximum up to t is clearly and 
hence is necessarily less than x m . The second term denotes 
the fact that in case of branching the maximum of each branch 
stays below x m : since the branches are independent, one gets 
a square. The third term corresponds to diffusion. By using 
(Ax) = and (Aa; 2 ) = 2Ddt and expanding Eq. [S] to the 
first order in dt and second order in Aa; we obtain 



dt 



dxl 



-y(Ro + l)Q + yR Q [9] 



for x m > 0, satisfying the boundary conditions Qt(0) = and 
Qt(oo) = 1, and the initial condition Qo(x m ) = 0(x m ), where 
<d is the Heaviside step function. Hence from Eq. [H] 



{L(t))=2n [1 - Qt(xm)]dx r , 



[101 



Equation [9] can be solved numerically for all t and all Ro, 
which allows subsequently computing (L(t)) in Eq. [To] (de- 
tails and figures are provided in the Supplementary Material) . 

The convex hull area. To compute the average area in Eq. M, 
we need to evaluate {x m (t)} as well as (i m ). Once the cumula- 
tive distribution Qt(x m ) is known, the second moment {x m (t)) 
can be directly computed by integration, namely, {x m (t)) = 
Jq 00 da; m 2a; m (l — Qt(x m )). To determine (t m ), we need to also 
compute the probability density pt(t m ) of the random vari- 
able tm- Unfortunately, writing down a closed equation for 



Pt(t m ) is hardly feasible. Instead, we first define Pt(x m ,tm) 
as the joint probability density that the maximum of the x 
component achieves the value x m at time t m , when the full 
process is observed up to time t. Then, we derive a backward 
evolution equation for Pt (x m , t m ) and then integrate out x m 
to derive the marginal density pt{t m ) — Pt{x m ,t m ) dx m . 
Following the same arguments as those used for Qt{x m ) yields 
a backward equation for Pt(x m ,t m ): 



P t+dt (x m ,t m ) = [1 - j(Rq + l)dt] (P t (x n 



Ax, t r , 



dt)) 



+2jR dtQt{x m )P t (x m ,t m - dt). [11] 

The first term at the right hand side represents the contribu- 
tion from diffusion. The second term represents the contribu- 
tion from branching: we require that one of them attains the 
maximum x m at the time t m — dt, whereas the other stays 
below Xm (Qt(xm) being the probability that this condition 
is satisfied). The factor 2 comes fro m the interchangeability 
of the particles. Developing Eq. 11 to leading order gives 



dt 



d 

dt m 



Pt 



D 



J? 

dx'r. 



■y(Ro + l) + 2jR Qt 



ft -[12] 



This equation describes the time evolution of Pt{x m ,t m ) in 
the region x m > and < t m < t. It starts from the ini- 
tial condition Po(x m ,tm) = 8(x m ) S(t m ) (since the process 
begins with a single infected with x component located at 
x = 0, it implies that at t = the maximum x m ~ and also 
t m =0). For any t > and x m > 0, we have the condition 
Pt(x m ,0) — 0. We need to also specify the boundary condi- 
tions at Xm = and x m — > 00, which read (i) Pt(oo,t m ) = 
(since for finite t the maximum is necessarily finite) and (ii) 
Pt{0,t m ) = 5(tm) qt(x m )\x m =o- The latter condition comes 
from the fact that, if x m = 0, this corresponds to the event 
that the x component of the entire process, starting at ini- 
tially, stays below in the time interval [0, t], which happens 
with probability qt(x m )\x m =o'- consequently, t m must neces- 
sarily be 0. Furthermore, by integrating Pt(xm,t m ) with re- 
spect to tm we recover the marginal density qt( x m ) - 

The numerical integration of the full Eq. [12] would be 
rather cumbersome. Fortunately, we do not need this. Since 
we are only interested in {t m ), it is convenient to introduce 



Tt(Xm) — / tmPt{Xm,tm)dt n 
JO 



[131 



from which t he a verage follows as (t m ) ~ J dx m Tt(xm)- Mul- 
tiplying Eq. 12 by t m and integrating by parts we get 



0_ 

dt 



T t - q t (xm) = 



dx\ 



D 1^2 + 27#oQt - 7 (#0 + 1) 



Tt, [14] 



with the initial condition Tq(x, 
ditions T t (0) = and T t (oo) = 
numerically, together with Eq. 
Supplementary Material), and 



= 0, a nd the boundary con- 
Eq. 



14 



can be integrated 
(details are provided in the 
e behavior of 



[151 



f- OO 

(A(t)) = 7T / dXm [2x m (l ~ Qt(x m )) ~ T t (x m )] 
JO 

as a function of time is illustrated in Fig. [2] 



The critical regime. We now focus on the critical regime 
Ro — 1. We begin with the average perimeter: for Ro = 1, 
Eq. [9] admits a stationary solution as t — > 00, which can 
be obtained by setting dQ/dt — and solving the resulting 
differential equation. In fact, this stationary solution was al- 
ready known in the context of the genetic propagation of a 
mutant allele [22] . Taking the derivative of this solution with 
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respect to x m , we get the stationary probability density of the 
maximum x m 



Qoc(Xm) — 9x 7Tl Q oo (Xm) — 



[161 



The average is (x m ) = J °° x m qoo{x m ) dx m = ^/6D/j, which 
yields then Eq. [I] for the average perimeter of the convex 
hull at late times. 

To compute the average area in Eq. [7], we need to 



also evaluate the second moment (x m (t)), which diverges as 
t — > ao, due to the power-law tail of the stationary probabil- 
ity density q ao (x m ) oc x^ 3 for large x m . Hence, we need to 
consider large but finite t. In this case, the time dependent 
probability density qt(x m ) displays a scaling form which can 
be conveniently written as 



q t {x m ) ~ qoa(x m )f 



[17] 



where f(z) is a rapidly decaying function with f(z - C 1) ~ 1, 
and Hz ^> 1) ~ 0. Using the scaling form of Eq. [17] and 
Eq. 9 one can derive a differential equation for f(z). But it 
turns out that we do not really need the solution of f(z). 

From Eq. |17| we see that the asymptotic power-law decay 
of qt{x m ) for large x m has a cut-off around x* m ~ yDt and 
f(z) is the cut-off function. The second moment at finite but 
large times t is given by (x m (t)) — J °° x m qt(x m ) dx m . Sub- 
stituting the scaling form and cutting off the integral over x m 
at x* m — cyt (where the constant c depends on the precise 
form of f(z)) we get, to leading order for large t, 



.(*)> 



(Xm) dx 

77 



(>D 



lnt 



[181 



Thus, interestingly the leading order result is universal, i.e, 
independent of the details of the cut-off function f(z) (the in- 
dependence is only in the subleading term). To complete the 
characterization of (A(t)) in Eq. [7], we still need to determine 
(t m ): in the Supplementary Material we explicitly determine 
the stationary solution Poo(x m ,tm) for Rq = 1. By following 
the same arguments as for (x m (t)), we show that 



{tm) 
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hit 



[191 



for large t, which leads again to a l ogarithmic divergence in 
time. Finally, substituting Eqs. 18 and [l9] in Eq. [7] gives 
the result announced in Eq. pi. 

A deeper understanding of the statistical properties of the 
process would demand knowing the full distribution Prob(L, t) 
and Prob(j4, t) of the perimeter and area. These seem rather 
hard to compute, but one can obtain the asymptotic tails 
of the distributions by resorting to scaling arguments. Fol- 
lowing the lines of Cauchy's formula (see the Supplementary 
Material), it is reasonable to assume that for each sample the 
perimeter scales as L(t) ~ x m (i). We have seen that the distri- 
bution of Xm(t) has a power-law tail for large t: q cxi (x m ) ~ x^ 
for large x m . Then, assuming the scaling L(t) ~ x m (t) and 
using Prob(L,£ — > 00) dL ~ qoo(x m ) dx m , it follows that at 
late times the perimeter distribution also has a power-law 
tail: Prob(L,i — > 00) ~ L~ A for large L. Similarly, using 
the Cauchy formula for the area, we can reasonably assume 
that for each sample A(t) ~ x m (t) in the scaling regime. Once 
again, using Prob(A, t — > 00) dA = q 00 (x m ) dx m , we find that 
the area distribution also converges, for large t, to a stationary 
distribution with a power-law tail: Prob(y4, t — > 00) ~ A~ 2 for 



large A. Moreover, the logarithmic divergence of the mean 
area calls for a precise ansatz on the tail of the area distribu- 
tion, namely, 



Prob(A,t) 



2AnD ,_ 2 , 
-»■ A h 



5o 



[201 



where the scaling function h(z) satisfies the conditions h(z <C 
1) = 1, and h(z 3> 1) — 0. It is not difficult to verify that this 
is the only scaling compatible with Eq. [2] . These two results 
are consistent with the fact that for each sample typically 
A(t) ~ L 2 (t) at late times in the scaling regime. Our scaling 
predictions are in agreement with our Monte Carlo simula- 
tions (see Fig. [2j|. The power-law behavior of Prob(j4, t) im- 
plies that the average area is not representative of the typical 
behavior of the epidemic area, which is actually dominate d by 
fluctuations and rare events, with likelihood given by Eq. 20 



The supercritical regime. When Ro > 1, it is convenient to 
rewrite Eq. [9] in terms of W(x m ,t) = 1 — Q(x m ,t): 



5" 



D 



-W + j{R - l)W--yR W 2 



[21] 



starting from the initial condition W(x m , 0) 



(see Fig. From Eq. 



10 



<£(*)) 



; for all x m > 
W(x m ,t) dx m is 
m , up to a factor 



2-/o 

just the area under the curve W(x m ,t) vs. x 
2-7T. As t — > 00, the system approaches a stationary state f or al l 
Ro > 1, which can be obtained by setting dt W — in Eq. pi] . 
For Ro > 1 the stationary solution W(x m , 00) approaches the 
constant 1 — 1/Ro exponentially fast as x m — > 00, namely, 
W(x m ,oo) — 1 + -fi^ 1 — > exp[— x m /£], with a characteristic 
length scale ^ = D/ r y{Ro — 1). However, for finite but large 
t, W(x m ,t) as a function of x m has a two-step form: it first 
decreases from 1 to its asymptotic stationary value 1 — 1 /Ro 
over the length scale £ , and then decreases rather sharply from 
1 — 1 /Ro to 0. The frontier between the stationary asymptotic 
value 1 — 1 /Ro (stable) and (unstable) moves forward with 
time at constant velocity, thus creating a travelling front at the 
right end, which separates the stationary value 1 — 1/Ro to the 
left of the front and to the right. This front advances with 
a constant velocity v* that can be estimated using the stan- 
dard velocity selection principle |27l 1281 129] . Near the front 
where the nonlinear term is negligible, the equation admits 
a travelling front solution: W(x m ,t) ~ exp[— A(a; m — vt)], 
with a one parameter family of possible velocities v(\) = 
D\ + "/(Rq — 1)/A, parametrized by A. This dispersion rela- 
tion v(\) has a minimum at A : 
v(\*) 



A* = ^tXRo - l)/D, where 
1). According to the standard ve- 
[271 EH [29], for a sufficiently sharp 



v = V(A ) — 2^Dj(R 
locity selection principle 
initial condition the system will choose this minimum velocity 
v* . The width of the front remains of ~ 0(1) at large t. Thus, 
due to this sharpness of the front, to leading order for large 
t one can approximate W(x m ,t) ~ (1 — l/Ro)Q(v*t — x m ) 
near the front. Hence, to leading order for large t one gets 
(x m (t)) ~ (1 - 1/R )v*t and (x m ) ~ (1 - 1/Ro) (v'tf. The 
former gives, from Eq. [fjj^the result announced in Eq. 
For the mean area in Eq. l7J, the term (x m ) ~ t 2 for large t 
dominates over (t m ) ~ t (which can be neglected), and we get 
the result announced in Eq. pj. 

Conclusions 

In this paper, we have developed a general procedure for as- 
sessing the time evolution of the convex hull associated to the 
outbreak of an epidemic. We find it extremely appealing that 
one can successfully use mathematical formulae (Cauchy's) 
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from two dimensional integral geometry to describe the spa- 
tial extent of an epidemic outbreak in relatively realistic sit- 
uations. Admittedly, there are many assumptions in this epi- 
demic model that are not quite realistic. For instance, we have 
ignored the fluctuations of the susceptible populations during 
the early stages of the epidemic: this hypotheses clearly breaks 
down at later times, when depletion effects begin to appear, 
due to the epidemic invading a thermodynamical fraction of 
the total population. In addition, we have assumed that the 
susceptibles are homogeneously distributed in space, which is 
not the case in reality. Nonetheless, it must be noticed that 
in practical applications whenever strong heterogeneities ap- 
pear, such as mountains, deserts or oceans, one can split the 
analysis of the evolving phenomena by conveniently resorting 
to several distinct convex hulls, one for each separate region. 
For analogous reasons, the convex hull approach would not be 
suitable to characterize birth-death processes with long range 
displacements, such as for instance branching Levy flights. 

The model discussed in this paper based on branching 
Brownian motion is amenable to exact results. More gener- 
ally, realistic models could be taken into account by resorting 
to cumbersome Monte Carlo simulations. The approach pro- 
posed in this paper paves the way for assessing the spatial 
dynamics of the epidemic by more conveniently solving two 
coupled nonlinear equations, under the assumption that the 
underlying process be rotationally invariant. 

We conclude with an additional remark. In our computa- 
tions of the mean perimeter and area, we have averaged over 
all realizations of the epidemics up to time t, including those 
which are already extinct at time t. It would also be interest- 
ing to consider averages only over the ensemble of epidemics 
that are still active at time t. In this case we expect differ- 
ent scaling laws for the mean perimeter and the mean area of 
the convex hull. In particular, in the critical case, we believe 
that the behavior would be much closer to that of a regular 
Brownian motion. 

Appendix: Supplementary Materials 

Cauchy's formula. The problem of determining the perime- 
ter and the area of the convex hull of any two dimensional 
stochastic process [x(t), y(r)] with < r < t can be mapped 
to that of computing the statistics of the maximum and the 
time of occurrence of the maximum of the one dimensional 
component process x(t) |23l 124] , This is achieved by resort- 
ing to a formula due to Cauchy, which applies to any closed 
convex curve C. 

A sketch of the method is illustrated in Fig. 5. Choose 
the coordinates system such that the origin is inside the curve 
C and take a given direction 9. For fixed 9, consider a stick 
perpendicular to this direction and imagine bringing the stick 
from infinity and stop upon first touching the curve C. At 
this point, the distance M{9) of the stick from the origin is 
called the support function in the direction 9. Intuitively, the 
support function measures how close can one get to the curve 
C in the direction 9, coming from infinity. Once the support 
function M{9) is known, then Cauchy's formulas [30] gi ve the 
perimeter L and the area A enclosed by C, namely 

L = ^M(6)dB 
A=\S^[M\9)-(M>{9)f}d9, [22] 

where M 1 (9) — dM/d9. For example, for a circle of 
radius R = r, M{9) = r, and one recovers the stan- 
dard formulae: L = 2ivr and A — ivr 2 . When C is 
the convex hull of associated with the process at time t, 
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we first need to compute its associated support function 
M{9). A crucial point is to realize that actually M(9) = 
maxo<i-<t [x(t) cos(f?) + y(r) sin((9)] |23l 124] . Furthermore, if 
the process is rotationally invariant any average is indepen- 
dent of the angle 9. Hence for the average perimeter we 
can simply set 9 = and write (L(t)) = 27r(M(0)}, where 
brackets denote the ensemble average over realizations. Sim- 
ilarly for the average area, (A(t)) = vr [(M 2 (0)} - (Af'(0) 2 }]. 
Clearly, M(0) = maxo< T <t[i(r)] is then the maximum of the 
one dimensional component process x(t) for r £ [0,t]. As- 
suming that x(t) takes its maximum value x(t m ) at time 
t = t m (see Fig. 4). Then, M(0) = x(t m ) = x m (t), and 
M'(0) = y(t m ) [31]. Now, by taking the average over Cauchy's 
formulas, and using isotropy, we simply have Eqs. [5] and [6] 
for the mean perimeter and the mean area of the convex hull C 
at time t. Note that this argument is very general and is appli- 
cable to any rotationally invariant two dimensional stochastic 
process. Since the branching Brownian motion with death is 
rotationally invariant we can use these formulae. 

Numerical methods. Numerical integration. Equations [9] and 
[14] have been integrated numerically by finite differences 
in the following way. Time has been discretized by setting 
t = ndt, and space by setting x = idx, where dt and dx are 
small constants. For the sake of simplicity, here we consider 
the case Ro = 1. We thus have 

Qn + l(i) = 

= Q n (i) +fdt[l ~Q n (i)] 2 + 
Dj^[Q n (i + l)-2Q n {i) + Q n (i-l)] [23] 

and 

Tn + l(l) = 

= T n (i) + 2jdtT n (i) [Q n (i) - 1] + 

D Td3p [T " (l + ^ ~ 2T ' iW + Tn< - 1 ~ 1)] + 

£[T n (i)-T„(i-l)]. [24] 

As for the initial conditions, Qo(0) = and Qo{i > 0) = 1, 
and To(i) = Vi. The boundary conditions at the origin 
are Q n (0) = and T n (0) = 0. In order to implement the 
boundary condition at infinity, we impose Qn(i ma x) = 1 and 
Tii(jmax) = Vn, where the large value i max is chosen so that 
rn(imax)— Tn(*max — 1) < 10~ 7 . We have verified that numeri- 
cal results do not change when passing to the tighter condition 
T n {i max) T n (imax - 1) < 10" M . 

Once Q n (i) and T n (i) are known, we use Eqs. [10] and [15] 
to determine the average perimeter and area, respectively. 

Monte Carlo simulations. The results of numerical in- 
tegrations have been confirmed by running extensive Monte 
Carlo simulations. Branching Brownian motion with death 
has been simulated by discretizing time with a small dt: in 
each interval dt, with probability bdt the walker branches and 
the current walker coordinates are copied to create a new ini- 
tial point, which is then stored for being simulated in the 
next dt; with probability ydt the walker dies and is removed; 
with probability 1 — (b + 7) dt the walker diffuses: the x and 
y displacements are sampled from Gaussian densities of zero 
mean and standard deviation \/2Ddt and the particle posi- 
tion is updated. The positions of all the random walkers are 
recorded as a function of time and the corresponding convex 
hull is then computed by resorting to the algorithm proposed 
in [32]. 

PNAS I Issue Date | Volume | Issue Number | 5 



Perimeter statistics. In order the complete the analysis 
of the convex hull statistics, in Fig [6] and Fig [7] we show the 
results for the perimeter. 

Analysis of t m .In the critical case Ro = 1, the stationary 
joint probability density Pao(x m ,t m ) satisfies (upon setting 
dPt/dt = in Eq. [12]) 



d 
dt Wl 



D 



2l 



Poo \psm i t"n 



Poc \pCrn i in 



[251 



For any x m > 0, we have the con dition Poo(a;m,0) = 0. The 
boundary conditions for Eq. |25] > are Pao{x m — > oo,t m ) = 
and P x (0,t m ) = qoo {0)8{t m ) = 2y/-y/(6D)6(t m ). We first 
take the Laplace transform of (|25[), namely, 



Poo y^m 5 



This gives for all x m > 



s dx 2 K 



p (• i \dt 



[261 




[27] 

where we have used the condition Poo(x m ,0) = for any 
x m > 0. This second order differential equation satisfies 
two boundary conditions: Poo(oo,s) = and Poo(0, s) = 
2-^/77(613). The latter condition is obtained by Laplace trans- 
forming Poo(0, t 

77 



2^j/(6D) 5(t m ). By setting 



PootCs) = 2 v /7/(6£>). 
finally get 



By reverting to the variable x m , we 




Now, we are interested in determining the Laplace trans- 
form of the marginal density Poo(s) = J °° e~ stm Pcx,(t m ) dt m 
where p^^m) = J °° Poo(x m ,t m ) dx m . Taking Laplace trans- 
form of this last relation with respect to t m gives Poo(s) = 
J °° Poo{x m , s) dx m . Once we know Poo(s), we can invert it to 
obtain p x (tm)- Since we are interested only in the asymptotic 
tail of Poo(im), it suffices to investigate the small s behavior 
of Poo(s). Integrating Eq. (31 1 over x m and taking the s — > 



limit, we obtain after some straightforward algebra 



Poo(s) = 1 + — s ln(s) + 

07 

We further note that 



2 / \7 d ~. , \ 3 

tin Poc{tm) dtm = , Poc(S) — — , 

ds^ 07s 



[321 



[331 



which can then be inverted to give the following asymptotic 
behavior for large t m 



Poo (tm ) 



[341 



Analogously as for {x m ), the moment (t m ) — > 00, due to the 
power-law tail Poo(i m ) oc t m 2 - Hence we need to compute (t m ) 
for large but finite t: in this case, the time-dependent solution 
displays a scaling behavior 



+ X T , 



we rewrite the equation as 
^P - 

OZ l 



^Pao = 0. 



[281 



[291 



Upon making the transformation Poo(z) = \[z F(z), the func- 
tion F(z) then satisfies the Bessel differential equation 



d 2 „, , 1 d . 
dzT^ + -zdz F ^ 



1 + 



49 

4? 



F{z) = 0. 



[301 



The general solution of this differential equation can be ex- 
pressed as a linear combination of two independent solutions: 
F(z) — AI 7 / 2 (z)+B K 7 / 2 (z) where I v (z) and K v (z) are mod- 
ified Bessel functions. Since, I v (z) ~ e z for large z, it is 
clear that to satisfy the boundary condition Poo(oo,s) = 
(which mean F(z — > 00) = 0), we need to choose A = 0. 
Hence we are left with F(z) = BK 7 / 2 (z), where the con- 
stant B is determined from the second boundary condition 



Pt (tm) — Poo(tm) g 



[351 



where the scaling function g(z) satisfies the conditions g(z <C 
1) ~ 1 and g(z > 1) = 0. Much like in Eq. [17] for the 
marginal density q t (x m ), we have a power-law tail of pt(t m ) 
for large t m that has a cut-off at a scale t* m ~ t, and g(z) is 
the cut-off function. As in the case of x m , we do not need the 
precise form of g(z) to compute the leading term of the first 
moment (t m ) = J Pt(t m ) t m dt m for large t. Cutting off the 
integral at t m = c\t (where C\ depends on the precise form of 
g{z)) and performing the integration gives 



(tm) 



tm Poo(^m) dtn 



57 ' 



lnt, 



[361 



which is precisely the result announced in Eq. [19] 
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Fig. 1. The snapshots of the trajectories of an assembly of infected individuals at the epidemics outbreak at three different times (schematic), starting from a single infected 
at the origin O at time t = 0. Individuals that are still infected at a given time t are displayed as red dots, while those already recovered are shown as black dots. The convex 
hull enclosing the trajectories (shown as a dashed line) is a measure of geographical area covered by the spreading epidemic. As the epidemic grows in space, the associated 
convex hull also grows in time. 
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Fig. 2. Left. The average area (A(t)) of the convex hull as a function of the observation time. For the parameter values, we have chosen D = 1/2 and = Ro'f = 0.01. 
We considered five different values of Rq. We have obtained these results by two different methods: (i) via the numerical integration of Eqs. |9j and [14] and using Eq. [15] . 
These results are displayed as solid lines, (ii) by Monte Carlo simulations of the two-dimensional branching Brownian motion with death with the same parameters, averaged 
over 10 samples. Monte Carlo are displayed as symbols. The dashed 
the numerical simulations are provided in the Supplementary Material. 



lines represent the asymptotic limits as given in Eq. Jij for the critical case Rq = 1. Further details of 
Right. Distribution of the area of the convex hull for the critical case Rq = 1, with "f = 0.01 and 



D = 1/2, as obtained by Monte Carlo simulations with 2 ■ 10 realizations. The dashed line corresponds to the power- law (247rD/57)A 2 as predicted by Eq bol. 
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Fig. 3. Left. The time behavior of the average area in the supercritical regime for different values of Rq > 1. Dashed lines represent the asymp toti c scaling as in Eq. 
[4] . The red curve corresponds to the critical regime. Right. The behavior of Wt{x) = 1 — Qt{x) for Rq = 1.5 at different times, as in Eq. [2l] , When t — > 00, 
Wt (x) — > 1 — Rq , and for large but finite times the travelling front behavior is clearly visible. The inset displays the exponential convergence of Wt (x) to the asymptotic 
limit. The dashed line represents £ = y/D /y^Ro — 1). 
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Fig. 4. Left. A branching random walk composed of five individuals. At time t = 0, a single infected is at the origin O, and starts diffusing (blue line). At later times, this 
individual branches and gives rise to other infected individuals. Among these, the red path reaches the maximum x m along the X component up to the final time t. Infected 
individuals at a given time t are displayed as red dots, whereas recovered as black dots. Center. The displacement along the X direction as a function of time. The red path 
reaches the global maximum x m at time t m . Right. The displacement along the y direction as a function of time. When the red path reaches the global maximum £ m 
at time t mi its y coordinate attains the value yf£ m ). A crucial observation is that the y component of the trajectory connecting O to the red path is a regular Brownian 
motion. This is not the case for the X component, which is constrained to reach the global maximum of the branching process. 




Fig. 5. Cauchy's construction of the two-dimensional convex hull, with support function 
M{9) representing the distance along the direction 6. 
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Fig. 6. Left. The average perimeter (L(t)) of the convex hull as a function of the observation time. For the parameter values, we have chosen D = 1/2 and 
b = Hoy = 0.01. We considered five different values of Rq. We have obtained these results by two different methods: (i) via the numerical integration of Eq. [9] and using 
Eq. [10] (with the choices dt = 0.003125 and dx = 0.1768). These results are displayed as solid lines, (ii) by Monte Carlo simulations of the two-dimensional branching 
Brownian motion with death with the same parameters and with the choice of the Monte Carlo time step dt = 0.25 with the results averaged over 10 samples. Monte 
Carlo are displayed as symbols. The dashed lines represent the asymptotic limits as given in Eq. [1] for the critical case Rq = 1. Right. Distribution of the perimeter of the 
convex hull for the critical case Rq = 1, with 7 = 0.01 and D = 1/2, as obtained by Monte Carlo simulations with time step dt = 1 and t = 4 ■ 10 . The number of 
realizations is 2 ■ 10^. The dashed line of the left panel corresponds to the power-law L~ (up to an arbitrary prefactor). 
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